function [time, S] = GPM_Cyclic(C, S0, m)

[N, ~] = size(S0);
d = 2;
n = N/d;

iter_max = 100;
tol = 1e-8;

time_start = tic;

% Initialization
S = S0;
iter = 1;
val_last = inf;
Y = C*S;
val = trace(S'*Y);

while ( iter <= iter_max && abs(val - val_last)/N > tol ) 
%     tic
    for i = 1:n           
        S( (i-1)*d+1:i*d , : ) = Proj_Cyclic(Y((i-1)*d+1:i*d,:),m);
    end
    fprintf('GPM iter = %d\n', iter);
    Y = C*S;
    val_last = val;
    val = trace(S'*Y);
    iter = iter + 1;
end

% S = S*S(1:d,:)';

time = toc(time_start);